Influence of Spatial Correlations on the Lasing Threshold of Random Lasers 
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The lasing threshold of a random laser is computed numerically from a generic model. It is shown 
that spatial correlations of the disorder in the medium (i. e., dielectric constant) lead to an increase 
of the decay rates of the eigenmodes and of the lasing threshold. This is in conflict with predictions 
that such correlations should lower the threshold. While all results are derived for photonic systems, 
the computed decay rate distributions also apply to electronic systems. 



In the theory of disordered media two important 
regimes, diffusive and localised, are distinguished [I). In 
the diffusive regime (for weak or moderate disorder), 
eigenstates are extended and efficient transport is pos- 
sible. In the localised regime (for strong disorder) , eigen- 
states become localised and transport is strongly inhib- 
ited. Many experimental findings for random lasers are 
more consistent with the assumption of a lasing mode 
that is localised while a direct experimental analysis of 
the sample shows that it is in the diffusive regime. 

To determine whether a sample is in the localised or 
in the diffusive regime, a transport property is measured. 
The most efficient way to achieve this is to check for the 
rounding of the backscattering cone Q. Such a round- 
ing is not reported from experiments 0, Q . Transport 
is, however, dominated by extend eigenstates, and the 
simultaneous existence of a few localised eigenstates in 
a sample in the diffusive regime, i.e., that is on average 
diffusive, would not be noticed pg. Such localised modes 
have recently been detected experimentally in a diffusive 
sample |fg. 

The important question is to explain under which con- 
ditions such localised eigenstates can exist in a diffusive 
sample. [These states have been termed anomalously lo- 
calised states (ALS) or prelocalised states. For a recent 
review see Ref. 7>] One-dimensional disordered systems 
are always in the localised regime, i.e., they can never 
show diffusive behaviour. Theoretical studies on such 
systems thus cannot give information on the interplay 
between extended and localised modes. The situation is 
different in two- and three-dimensional samples. Two- 
dimensional samples shorter than the localisation length 
behave similarly to three-dimensional samples, and one 
is allowed to replace three-dimensional systems by their 
computationally cheaper two-dimensional counterparts. 

The computational cost of treating a two-dimensional 
sample is significantly higher than for a one-dimensional 
sample, and only few studies have been published [8j- 
Reference models the scatterers in the disordered me- 
dia as dipoles, Ref. |9| studies circular particles using 
the finite-difference time-domain (FDTD) method. Both 
publications do not state explicitly whether their sam- 
ples are in the diffusive or in the localised regime but the 
parameters given strongly suggest that the samples are 
in the localised regime. 

The only publication so far on the interplay between 
diffusive and localised eigenstates inside a diffusive sam- 
ple seems to be Ref. ^3- There it was estimated that lo- 



calised states become exponentially more frequent when 
the disorder inside the sample is spatially correlated. Mo- 
tivated by the picture of photons travelling in closed loop 
inside a ring-shape structure [3| , they study a ring-shaped 
area of higher dielectric constant. This is a very special 
situation, and it is not obvious how characteristic such a 
special situation is for the entire behaviour. (It should be 
noted that the opposite effect, namely that in a localised 
sample a few modes become extended when spatial corre- 
lations in the disorder are introduced, is well-understood, 
see e. g. Ref. ITU) 

In this paper we will study this problem from a more 
generic approach. The lasing threshold of a sample is 
determined by the decay rates of the eigenstates of the 
system since the loss (=decay) of photons in the mode 
has to be compensated by pumping if the sample is to 
start lasing action. Following the approach of Ref. ^3 
we numerically compute the decay rate distribution of 
a two-dimensional sample on a suitable grid. (Earlier 
work on the lasing threshold of chaotic cavities [Ij] can- 
not be applied since by construction all eigenstates are 
extended |l|.) We improve on previous work by including 
spatial correlations. 

We use the Anderson Hamiltonian which describes the 
motion of an uncharged particle in a spatially varying 
potential. The results can directly be applied also to 
photonic systems since the Helmholtz equation with a 
spatially varying dielectric constant has the same form. 
The sample is discretised with lattice spacing A, where 
for electronic systems A = Tt/k-p (kp is the wave vector 
at the Fermi level) and for photonic systems A = 2\/ir 
(A is the wave length of the light). This is a natural 
choice in that there is then one propagating mode per 
transversal lattice point, and the width of the sample is 
best measured in terms of the number N of propagating 
modes. 

Transport is modelled by nearest-neighbour hopping 
with rate 1. (The results are easily adapted to arbitrary 
speed c of transport.) With a spatially varying potential 
P{x, y) the Hamiltonian for a sample of length L = LA 
becomes ^2 



'H(x, y ).(x\y') = S xx 'S yyl [P(x,y) - i(Si y + S L y)} 

+ 8yyl (<5a;+l,X' + 5'X-1,X') + 8 XjX ' {5y+l,y> + Sy-l.y') 



(1) 



with x = 1, . . . ,L and y = 1, . . . , N . The imaginary 
part of Ti models coupling of the sample to the outside 
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where we assume that we operate at the centre of the 
conduction band. 

The spatial correlations are assumed to fall of exponen- 
tially such that P(x, y) takes on random values, normal- 
distributed with zero mean and correlator 

(P(f)P(r')) = u, 2 exp(-^^) . (2) 

w measures the strength of the disorder, and R c is 
the correlation radius. Since we need to generate a 
large number of mutually correlated random numbers, 
a Fourier based method has to be employed 

The eigenvalues of the matrix Ti correspond to the 
(quasi-)eigenmodes of the system. Their real part u> gives 
the energy (or, for photonic systems, the frequency) of 
the mode, and their imaginary part 7 the decay rate |l5| . 
We thus have an eigenvalue problem of a non-Hermitian 
complex symmetric matrix but an eigensolver specifically 
adopted to this structure exists [12j. Even with this ef- 
ficient eigensolver, this is still a numerically expensive 
task, and it is impossible to analyse so many samples 
that there would be no more noise in the results. 

While the model is described in terms of the disorder 
strength w, contact with experiments or analytical theo- 
ries is best made by introduction of the mean-free path 
/. It can be computed from the length-dependence of the 
transmission probability T through the sample. In the 
diffusive regime, I < L <C Nl, it is given by [l| 



The transmission probability has been computed using 
the method of recursive Green's functions for vari- 
able disorder strength w and correlation length R c . We 
determined the mean-free path by fitting the numerically 
computed T(L) to this functional form self-consistently 
in the interval [/; 101]. (Picking some other interval, e. g. 
[0; I] , changed the result only by about 1 %.) The rescaled 
results are depicted in Fig. ^ f° r both N — 51 and 
TV = 81, i.e., for samples of different width. As the figure 
shows, both sets of curves are almost identical, thereby 
demonstrating that we are operating in the wide-sample 
regime. The mean-free path increases significantly as R c 
increases. This is immediately obvious since with increas- 
ing R c the potential changes less within a given distance; 
hence, there is less scattering. 

We would like to point out two "curiosities" . The nu- 
merical data suggest that the mean-free path I factorises 
as l(w,R c ) = fi(w)f2(Rc)- We did not manage to find 
an explanation for this observation. Furthermore, the 
mean-free path seems to scale as I ex iu — 1 ' 71 where 1.71 
is a numerical parameter. For uncorrelated random or- 
der that is uniformly distributed in the interval [—w; w], a 
scaling I oc to -1-5 was found numerically [l^. An analyt- 
ical theory is available only for one-dimensional systems 
in the limit w — ► where I oc 1/w 2 is found , so that 
a universal scaling for finite w might not exist at all. 

The increase of I with increasing R c poses a problem 
for a systematic study of the effects of correlations: One 
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FIG. 1: Numerically computed rescaled mean-free path I, de- 
pending on the disorder strength w and the correlation radius 
Rc (in units of the lattice spacing A). The solid lines are for 
samples of width N = 51, the dashed lines for samples of 
width N — 81. Samples were computed with w in steps of 
0.1 and R c in steps of 0.5 A (plus the value 7? c = 0.2 A). By 
rescaling I — > Iw 1 we can demonstrate the apparent scaling 
I tx w~ 1,71 and the factorisation l(w,R c ) = fi(w)f2(R c )- 



has to decide whether to compare samples with identical 
I (and thus variable w) or samples with identical w (and 
thus variable I). The final results must depend (apart 
from trivial prefactors) only on the ratios L/l and R c /l — 
not on any of those quantities separately. This decision is 
thus "only" one of numerical efficiency and minimisation 
of finite-size effects. 

For most of our simulations, we have decided to keep I 
constant at I — 12.5 A. For each value of R c , the needed 
value for w was determined by interpolation of the nu- 
merical data presented in Fig. ^ The choice of constant 
I offers the advantage that, even if R c is changed, sam- 
ples with identical "physical" length L/l occupy the same 
number of lattice points, and thus need the same amount 
of numerical work. (For constant "physical" length L/l, 
the needed computing time scales as 0(l 2 ). With con- 
stant w, this would impose severe restrictions on the 
range of R c that could be treated.) 

We have computed the decay rates for samples of width 
N = 50 for length L/l = 1,2,3,4,5,6,9,12,15,18, and 
correlation radius R c /A = 0.0, 0.2, 0.5, 1.0, 1.5, 7.5. 
For each set of parameters, approximately 2000 samples 
were generated. The maximum value of L is limited be- 
cause we are interested in the diffusive regime, hence L 
has to be sufficiently smaller than the localisation length 
Lioc = (N +1)1/2. We did not consider larger values of R c 
than 7.5 A since the sample should be much wider than 
the characteristic length scale of the disorder. Otherwise 
the sample would effectively become one-dimensional. 

To check the results, we have computed the decay rate 
distribution also for N — 80 for a few selected values of 
L/l and R c . To complement the other simulations, we 
have kept w constant. As explained above, this implies 
that we could only include R c < 2A. 

Following the approach introduced in Ref. ^3 f° r sam- 
ples in the diffusive regime, we fit the numerically com- 
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FIG. 2: Characteristic decay rate 70 as a function of sample 
length L and correlation radius R c (for samples of width N = 
50). The dashed line is for the control simulations with N = 
80 and L = 25 1. 



puted decay rate distribution to the functional form 

M7. 



1-Q(M + 1, 



7o 



(4) 



where the fitting parameters M and 70 depend on N, 
L and R c , and Q(a,x) = T(a,x)/T(a) is the regularised 
Gamma function . 

All numerically computed histograms fit well to the 
form (0J. The dependence of P("f) onto M is only weak 
for M ^> 1, making a precise determination of M diffi- 
cult. Within this error limit, we did not find a significant 
dependence of M on R c , and M is approximately given 
by the R c = result M = N/[1 + L/(6l)\ 

The fitting parameter 70, marking the typical value of 
the decay rates, can be determined to much better preci- 
sion. 70 is much more important for the lasing threshold 
than M, so the limited precision of M does not pose a 
problem. The determined 70 is shown in Fig. [5] 

For R c = 0, jqL 2 seems to approach a constant value 
as L is increased. This value is about 20 % larger than 
the value l/(2cZ) found numerically for equi-distributed 
disorder in the interval [— w;w] [12|. Trying to approach 
the limit L — > 00 numerically is not possible since then 
the sample would become localised. 

The important conclusion from Fig- E] is that for sam- 
ples of arbitrary length, the introduction of correlations 
in the disorder leads to an increase of the decay rates. 
This increase is quick as R c is increased starting from 0, 
and becomes slower for large R c . The same behaviour is 
seen in the control simulations with TV = 80 and fixed w 
(and thus variable I). 

Until now, all results are valid for both electronic and 
photonic systems. Now we will specialise to random 
lasers. The light inside a random laser is amplified by 
a laser dye. This dye is able to amplify light within a 
certain range of frequencies, so only K 3> 1 eigenmodes 
out of all eigenmodes of the system are amplified. This 
number varies only slightly between different realisations 
of the same ensemble due to an effect known as spectral 
rigidity 1] . The lasing threshold is given by the smallest 
decay rate out of the K modes within the amplification 



window . This is immediately obvious since the lasing 
threshold is passed when photons are created faster than 
they can decay (=escape from the sample). 

There are two different approaches to computing the 
lasing threshold of a random laser. The direct approach is 
to compute the eigenmodes of a certain number of realisa- 
tions of the disordered systems, then for each realisation 
to determine the smallest decay rate inside the amplifica- 
tion window, and finally collect statistics for those values. 
Since this process yields only a single datum per sample, 
a very large number of realisations needs to be computed 
to arrive at data of sufficient quality. The average lasing 
threshold determined in this way is depicted in Fig. as 
dashed line. 

Frequently more efficient is the second approach where 
one starts with the computation of the distribution P(-f) 
of the individual decay rates. The intermediary result is 
either a numerical histogram, or, by fitting the histogram 
to an analytical form, a distribution function that can be 
evaluated directly for arbitrary argument. We adopt the 
latter and use Eq. together with the values of M and 
70 computed by fitting. 

The distribution Pi (71) of the lasing threshold is the 
distribution of the smallest value out of K values, each 
distributed according to P(j). This assumes that the 
decay rates of different modes are uncorrelated. For K 
I this seems logical but to our best knowledge no explicit 
check of this assumption has been published so far. As a 
side-effect of our computations, we will fill this gap. 

Pi (71) is difficult to evaluate numerically for K ^> I 
since it is sharply peaked. The position 7 m of the maxi- 
mum of Pi is immediately seen to be given by 
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(^-l)[P(7 m )] 5 



(5) 

Since p is that sharply peaked, 7 m already contains 
all the relevant information, and nothing relevant is lost 
by not computing the entire distribution. The lasing 
threshold computed from Eq. (J5J, after inserting the fit- 
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FIG. 3: Comparison of the average of the lasing threshold, 
directly computed from the numerical data (solid line), and 
the most likely lasing threshold computed from the distribu- 
tion of the individual decay rates (dashed line). Both lines 
have been computed from the same samples, explaining the 
correlation of the noise in the two sets of lines. 
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ting parameters M and 70 computed from the numerical 
histograms into Eq. Q , is shown as solid line in Fig. [3] 

From the figure, two important conclusions can be 
drawn. First, the lasing threshold computed via the two 
separate methods agrees well. (The "noise" of the two 
sets of curves is correlated since the same raw data was 
used as input for both methods.) This means that the 
decay rates of different modes indeed are uncorrelated. 
Furthermore, also fitting the numerical histogram to the 
form Q is a valid procedure. 

The second conclusion — the heart of this paper — is 
that introducing spatial correlations into the disorder of 
a random laser increases the lasing threshold, in contra- 
diction to predictions [T^j ■ 

In this paper, we have thus arrived at two related — 
but not identical — results. We have shown that the 
decay rates increase if spatial correlations of the disorder 
are introduced. The computed decay rate distribution 
possesses the same form, just with different parameter, 
as earlier observed for diffusive samples with uncorrelated 
disorder 0|. This first result means that the "typical" 
eigenstates become more lossy. 

Our second result is that also the lasing threshold in- 
creases. This means that also the "special" eigenstates 
with lower-than-average loss, which are selected by mode 
competition to become the lasing modes, become more 
lossy. Even though we did not directly compute the spa- 
tial extend of the eigenstates, this still clearly demon- 
strates that no (pre-)localised eigenstates are formed by 
the introduction of spatial disorder. We thus fail to 
observe the prediction that such states should be cre- 
ated yJJ. 



There are several explanations for the difference be- 
tween our results and Ref. 0. One explanation is that 
a single ring-shaped area of increased dielectric constant 
does lead to the formation of a localised state, as sug- 
gested by the authors, but the influence of the disorder 
around that ring-shaped area significantly reduces this 
effect. Another, equally likely, explanation is that in our 
simulations we are only able to treat samples of finite 
size, with a finite number of eigenstates. The creation of 
a localised state may be an event that is so rare the we 
fail to see such an event occur in our finite-size simula- 
tions. On the other hand, the typical length scale is given 
by the area per lasing mode, measured experimentally to 
be a few 10 /1m 2 in two-dimensional ZnO films |19|. and 
our samples are larger than this. 

To give a definite answer on whether spatial correla- 
tions can explain the formation of localised states, more 
numerical studies are needed, preferably using different 
methods. Specialised but numerically efficient models H 
cannot incorporate spatial correlations of the dielectric 
constant. Two-dimensional FDTD simulations have al- 
ready been used to describe random lasers |E| . They need 
to make only minimal assumptions and they can be ex- 
tended to include arbitrary spatial correlations. FDTD 
simulations thus might be a good candidate but diffu- 
sive samples need to be larger than the localised samples 
studied so far. Given that FDTD is computationally very 
expensive, it is not obvious to us whether this would still 
be numerically feasible. 
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